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The effective increase of the critical density associated with the interaction of relativistically 
intense laser pulses with overcritical plasmas, known as self-induced transparency, is revisited for 
the case of circular polarization. A comparison of particle-in-cell simulations to the predictions of 
a relativistic cold-fluid model for the transparency threshold demonstrates that kinetic effects, such 
as electron heating, can lead to a substantial increase of the effective critical density compared to 
cold-fluid theory. These results are interpreted by a study of separatrices in the single-electron 
phase space corresponding to dynamics in the stationary fields predicted by the cold-fluid model. 
It is shown that perturbations due to electron heating exceeding a certain finite threshold can force 
electrons to escape into the vacuum, leading to laser pulse propagation. The modification of the 
transparency threshold is linked to the temporal pulse profile, through its effect on electron heating. 
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I. INTRODUCTION 

The optical properties of a plasma under the action 
of a relativistically intense laser pulse (intensity I > 
10 18 Wcm~ 2 for 1 um wavelength) are profoundly af- 
fected by nonlinearities in the corresponding laser-plasma 
interaction. In particular, the question of whether a pulse 
with the carrier frequency ujt, propagates in a plasma of 
electron density no can no longer be answered solely in 
terms of the critical density, 

n c = e n m e u> 2 L /e 2 , (1) 

where m e is the electron rest mass, — e is the electron 
charge, and eo is the permittivity of free space. By def- 
inition, a relativistically intense pulse accelerates elec- 
trons from rest to relativistic momenta within an optical 
cycle and, thus, the electron mass in Eq. (1) has to be 
corrected by the relativistic factor 7 = y/l + p 2 /m 2 c 2 , 
where p is the electron momentum. For a purely trans- 
verse wave propagating through a cold, homogeneous 
plasma, this relativistic factor can be related, by the 
conservation of canonical momentum, to the normalized 
amplitude of the wave vector potential ao = eA /(m e c), 
7 ~ \J\ + a 2 l /2 [1]. Therefore, one is forced to introduce 
an intensity-dependent effective critical density [2, 3] 

nf = <fi^n e . (2) 

According to Eq. (2), a relativistically intense laser 
pulse (ao > 1) can propagate through a nominally over- 
dense plasma, with electron density n c < n e < n° ff , 
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a phenomenon known as relativistic self-induced trans- 
parency (RSIT). Apart from its role as a fundamental 
process in laser-plasma interaction, RSIT is also inter- 
esting for applications, as it often determines the regime 
of efficient laser-target interaction. In the context of ion 
acceleration, for instance, RSIT can prevent efficient ion 
radiation-pressure-acceleration from thin targets [4-8] or 
laser-driven hole-boring in thicker ones [9, 10]. On the 
other hand, RSIT may enhance electron heating in the 
break-out afterburner acceleration mechanism, thus al- 
lowing for higher ion energies [11-13]. 

In this paper, we investigate RSIT in the case of a 
circularly polarized (CP) laser pulse with finite rise (or 
ramp-up) time r r and infinite duration, normally inci- 
dent onto a semi-infinite plasma with a constant density 
uq > n c , and a sharp interface with the vacuum. This 
configuration is of particular interest for ultrahigh con- 
trast laser interaction with thick targets. Unfortunately, 
the simple relation (2), derived assuming a purely trans- 
verse plane-wave and a homogeneous plasma of infinite 
extent, does not apply to this setting. The main reason 
for this is that the effect of the ponderomotive force (asso- 
ciated here with inhomogeneities along the propagation 
direction) becomes dominant and leads to a significant 
modification of RSIT threshold. Since the 1970's, sev- 
eral analytical studies, mostly within the framework of 
relativistic, cold-fluid theory [2], have been undertaken 
to investigate strong electromagnetic wave propagation 
through inhomogeneous plasmas [14-16], culminating in 
a derivation of a modified RSIT threshold which incor- 
porates boundary conditions at the plasma-vacuum in- 
terface [17, 18]. In order to establish contact with this 
line of previous work and to focus on the key physical 
mechanisms, we will restrict attention to immobile ions 
and one-dimensional geometry. 

Based on the assumptions stated above, the relativis- 
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FIG. 1. (color online) (a) Effective critical density as a function of the laser field amplitude ao as predicted by the simple relation 
(2) (dashed, black line). Threshold density n t h(ao) below which, according to the cold-fluid theory (c/. Sec. II), no standing 
wave solutions exist (solid, blue line). RSIT threshold as extracted from our PIC simulations (c/. Sec. IV) with two different 
pulse rise times: 0.25ti (error bars) and 4tl (triangular error bars), where tl = 2-k/oji is the laser period, (b) Schematic 
representation of the stationary solution predicted by the cold-fluid theory for the case of total reflection [regions (A) and (B) 
of panel (a)]. Shown are the electric field E x (x), vector potential of the standing wave |a(x)|, and ion (electron) density rn(x) 
\n B {x)\, see Sec. II for details, (c) Schematic representation of a typical case of pulse propagation in PIC simulations [for RSIT 
in region (B) or (C) of panel (a)]. Arrows indicate the direction of electron motion. See Sec. IV for numerical results. 



tic cold-fluid model predicts total reflection of the inci- 
dent pulse, if a certain density threshold n t h(ao) is ex- 
ceeded [17, 18] [see solid blue line in Fig. 1(a)]. The 
geometry of the stationary state predicted for uq > 
Hth(^o) is illustrated in Fig. 1(b): the ponderomotive 
force pushes the electrons deeper into the plasma, cre- 
ating a charge separation layer (CSL) and an electron 
density spike [henceforth referred to as compressed elec- 
tron layer (CEL)] at the edge of the plasma. Electrons 
in the CEL experience a strong electrostatic field (due 
to charge separation) , which balances the ponderomotive 
force. The density in the CEL is typically much higher 
than n® ff and, thus, pulse propagation is inhibited and a 
standing wave is formed. 

For plasma densities uq < nth(oo)j such stationary so- 
lutions cease to exist, and one enters the regime of RSIT. 
particle-in-cell (PIC) simulations [18, 19], however, indi- 
cate that light propagation in this regime is quite differ- 
ent from the traveling- wave solutions discussed earlier [3] . 
Although a CEL is initially formed, electrons at its edge 
escape toward the vacuum, leading to force imbalance 
and allowing the ponderomotive force to push the CEL 
deeper into the target. The situation is more reminiscent 
of hole-boring [10, 20] (albeit with immobile ions) with 
a penetration front moving deeper into the plasma with 
a constant velocity Vf, and a Doppler-shifted reflected 
wave [Fig. 1(c)] (see also Refs. [21-24]). 

In this work we show, through PIC simulations, that 
in the presence of electron heating, induced by the pulse 
finite rise time, such a propagation mechanism can be 
activated even for densities no > n t h(ao); see Fig. 1(a). 
The crucial role is again played by electrons at the edge 



of the CEL escaping toward the vacuum. However, it has 
been recently shown that, in the total reflection regime, 
electrons at the edge of the CEL cannot be forced to es- 
cape into the vacuum by infinitesimal perturbations [19]. 
To interpret our results, we are thus led to study the 
response of electrons at the edge of the CEL to finite 
perturbations. Studying the dynamics of a test-electron 
in the stationary fields predicted by the cold-fluid model 
for the CSL and vacuum, we show that electron escape 
to the vacuum is controlled by separatrices in the single- 
electron phase space. Moreover, we demonstrate that 
the perturbation threshold for unbounded motion (elec- 
tron escape) predicted by our analytical considerations is 
comparable to the attainable electron momentum due to 
heating (in the CEL), observed in our PIC simulations 
at the threshold for RSIT. Finally, we study the effect 
of laser pulse rise time on electron heating and on the 
observed modification of the RSIT threshold. 



This paper is organized as follows. In Sec. II we revisit 
some results of the stationary cold-fluid theory that mo- 
tivate the present study. In Sec. Ill we analyze the single- 
electron phase space (for motion in vacuum and charge 
separation layer), by determining equilibrium solutions 
(Sec. IIIB), studying their linear stability (Sec. IIIC), 
and determining separatrices of bounded and unbounded 
motion (Sec. HID). In Sec. IV we present our PIC simu- 
lation results and relate them to the analytical results of 
Sec. III. Finally, we discuss our findings and present our 
conclusions in Sec. V. 
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II. REVIEW OF RELATIVISTIC COLD-FLUID 
THEORY FOR RSIT 

Throughout the paper, all quantities are normalized 
to (so-called) relativistic units. In particular, velocity, 
time, and distance are normalized to the speed of light 
c, inverse laser frequency uj^ 1 , and inverse vacuum wave 
number fc^ 1 = c/ljl, respectively. Electric charges and 
masses are normalized to e and m e , respectively, and 
densities are normalized to the critical density n c . Fi- 
nally, electric fields are normalized to the Compton field 
E c = m e cuj L /e. 



A. Stationary cold plasma model 

In this section, we revisit the one-dimensional station- 
ary model proposed independently by Cattani et al. [17] 
and Goloviznin and Schcp [18] to describe the reflection 
of an incident relativistic CP laser pulse by a nominally 
overdense plasma with constant electron density n$ > 1 
and a sharp interface with vacuum. Our presentation 
follows Ref. [17]. 

We consider an incident CP laser pulse propagating 
along the x-direction with the vector potential 



A L (t,x) 



V2 



[cos(t — x) y + sin(£ — x) z] , (3) 



where y and i denote the unit vectors forming an or- 
thonormal basis in the plane transverse to the laser prop- 
agation direction. The pulse is incident from vacuum 
(x < 0) onto a semi-infinite plasma (x > 0). In this 
work, as in Refs. [17, 18], we will neglect ion motion. 

As outlined in the Introduction, we will consider sta- 
tionary solutions expressing the balance of the pondero- 
motive and electrostatic forces, achieved once a CSL of 
sufficient thickness x b is created, see Fig. 1(b). Assum- 
ing total reflection of the laser pulse by the plasma, the 
balance of the radiation (~ a|) and electrostatic pres- 
sures [~ (no Xb) 2 provides a rough estimate for the 
thickness of the CSL, 



(4) 



The exact expression for Xb and the limits of applicability 
of Eq. (4) are discussed below; see Eq. (17). 

In the following, we will look for stationary solutions 
with vector potential of the form 

A(t, x) = a{x) [cos (t + 9/2) y + sin (t + 9/2) z] , (5) 

where 9 accounts for the phase jump of the reflected wave 

A R (t, x) = — % [cos(* + x + 9) y + sm(t + x + 9) z] 
v2 

at x = Xb and will be computed below. In what follows, 
we will refer to the spatial function a(x) as the "vector 
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a/7 




dx 


dx ' 




d 2 (f> 
dx 2 


= n e - 




d 2 a 
dx 2 


_ / n e 
V 7 





potential." Note that, in the absence of plasma, we have 
a(x) — V2ao cos (x + 9/2). 

Modeling electrons as a relativistic cold fluid, as in 
Ref. [17], we seek stationary solutions satisfying the sys- 
tem of equations 



(6) 
(7) 
(8) 



Here, 4>(x) is the electrostatic potential, n e (x) is the 
electron density and the Lorentz factor is written as 
'y(x) — \J\ + a?{x) through conservation of transverse 
canonical momentum. Equation (6) expresses the bal- 
ance between the electrostatic and ponderomotive forces 
inside the plasma. Hence, it holds only for x > Xb- Equa- 
tion (7) is simply Poisson equation and Eq. (8) is the 
propagation equation (in the Coulomb gauge) for the field 
prescribed by Eq. (5). 

To solve the system of Eqs. (6)-(8), one considers the 
CSL and the CEL separately. The electron density n e (x), 
electrostatic field E x (x) — —d<p/dx and vector poten- 
tial a(x) are obtained in each layer. Solutions are then 
matched at the electron front x = Xb to ensure continuity 
of a(x), of its first derivative da/dx and of E x {x). 



B. Charge separation layer, < x < Xb 

The electrostatic field in the CSL, < x < Xb, is easily 
found by integrating Poisson Eq. (7) with n e — (no 
electrons) and boundary condition E x (0) = (to match 
the electrostatic field at the vacuum), 



E x (x) 



d<t> 
dx 



no x . 



(9) 



Thus, the electrostatic field for < x < Xb increases 
linearly, up to its maximum value Eb = E x (xb) = n a Xb- 
For total reflection at x = Xb we can integrate Eq. (8) 
once to get 



2ai 



'b ' 



(10) 



where a& = a(xb) is the vector potential at the plasma 
boundary. Here, we write the amplitude of the standing 
wave arising from the combination of the incident and 
reflected waves A L and A R , respectively, as 



a(x) = V2 



ao sin 



\V2 



ao 



(x - x b ) 



(11) 



which implies that in Eq. (5) we have 9/2 = tt/2 — 
arcsin(a;,/v / 2 ao) — x b - At this point, there are two un- 
known quantities, Xb and ab, which will be determined 



self-consistently by considering the region x > Xf,- Note 
that we assume ab > 0, while from Eq. (11) we have 
a'{xb) < so that the ponderomotive force dj/dx — 
7 _1 ada/dx pushes electrons deeper into the plasma, 
thus balancing the electrostatic force. 



C. Compressed electron layer, x > Xb 

We now derive equations for the electron density, vec- 
tor potential and electrostatic field in the plasma, x > Xb- 
Combining Eqs. (6) and (7), one can rewrite the normal- 
ized electron density in the plasma as a function of the 
vector potential a(x) and its first two derivatives, 



n e (x) =n Q + 



1 



vT 



d?a 
dx 2 



1 



(12) 



Substituting Eq. (12) in Eq. (8), we obtain a differential 
equation for the vector potential only: 



cPa 
dx 2 



r^(£) 2 -( 1+ " 2 -»^K < 13 » 



In the case of total reflection, Eq. (13) describes the 
evanescent field in the overdense plasma, and has to be 
solved with boundary conditions a(x) — > and da/dx — > 
for x — > +oo [25]. Equation (13) admits a first integral, 



1 (da 
2(1 + a 2 ) \d~x~ 



- l - (2n Vl + a 2 - a 2 ) = -no, (14) 



which may be used to derive a solution that satisfies the 
required boundary conditions [15], 



a(x) 



(n - 1) cosh [(x - x )/X s ] 
n cosh [(x — x )/X s ] — (n — 1) 



(15) 



where \ s = (no — 1)~ 1//2 is the classical skin-depth, and 
Xq is determined by ensuring the continuity of the vector 
potential at x = Xb- 

With a(x) inside the plasma provided by Eq. (15), one 
obtains n e (x) from Eq. (12), while Eq. (6) provides the 
electrostatic field in this region, 



E x {x) 



dx 



(16) 



Equation (16) together with Eqs. (9) and (10) and the 
continuity of the electrostatic field at x = Xb gives an 
explicit expression for the position of the electron front, 



Xb 




(17) 



Finally, from Eqs. (10) and (14), one obtains: 



This equation defines, for a given incident laser field am- 
plitude ao and initial plasma density no, the maximum 
evanescent field ab in the plasma. Solutions ab of Eq. (18) 
should satisfy the additional condition 



2 al - al > , 



(19) 



which follows from Eq. (10). 

Note that, in the limit 1 <C ab <C ao, Eq. (17) allows 
us to recover the approximate result Eq. (4). On the 
other hand, from Eqs. (17)-(18) we find that in the limit 
ab ^ I (correspondingly a 2 , <C no) Xb — 1a\jn' { j' 2 . 



D. Threshold for RSIT 

For a given plasma density no, Eq. (18) admits a so- 
lution only when the maximum evanescent field ab satis- 
fies [17] 



2(n + al) < 3 n ^l + a 2 . 



(20) 



As shown in Ref. [18], for n < 3/2, solutions com- 
patible with Eq. (19) can only be found in the region 
ao < 2no(no — 1). Thus, in this case, the threshold inci- 
dent laser amplitude reads 



a th = 2n o( n o - 1) 



(21) 



For no > 3/2 condition (19) is always fulfilled and 
Eq. (20) defines the regime of total reflection. The 
threshold for RSIT corresponds to equality in Eq. (20). 
The maximum evanescent field at the threshold then 
reads 



n -n 



, 3 / 9 „ 



n + 1 



(22) 



The threshold incident laser field amplitude a t ^ above 
which RSIT occurs in a plasma with initial density no 
is obtained by substituting ab = a# from Eq. (22) in 
Eq. (18), 



x th 



n (1 



1 + a | - 1 - a%/2 



(23) 



Depending on the density range, Eq. (21), respectively 
Eqs. (22)-(23), define a threshold amplitude ath(^o)j 
above which RSIT occurs, for a given plasma density. 
Alternatively, for a given incident amplitude ao, we may 
read Eq. (21), respectively Eqs. (22)-(23), as defining 
an effective critical density n t h(ao) below which RSIT 
occurs. This is illustrated in Fig. 1(a). Equation (21) 
yields 
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while Eqs. (22)-(23) can be inverted analytically in the 
limit ri{) 1, yielding 



9^00-12 , nth>l. (25) 



Thus, the asymptotic behavior of n t h(ao) in the limit 

1/2 

ao 3> 1 is nth a o > a J 11110 ! 1 more restricting condition 
than Eq. (2), which for large ao becomes n° ff oc ao. 

As discussed in the Introduction, our PIC simulations 
indicate that for pulses with finite rise time, the transi- 
tion between total reflection and RSIT occurs within the 
limits set by Eqs. (25) and (2) and, moreover, depends 
on the pulse rise time. In order to explain this discrep- 
ancy, we will now study single electron dynamics in the 
stationary fields (in vacuum and CSL) calculated above. 




III. SINGLE ELECTRON DYNAMICS 

A. Equations of single electron motion 

The equations of motion for an electron in the region 
x < xi, (i.e. in the vacuum and CSL), in the case of total 
reflection, read 



Px 



PxH, 
dx 



EJx) 



(26) 
(27) 



where we have used conservation of transverse canonical 
momentum to write the electron 7 factor as 



j(x,p x ) = \f\ + a 2 (x) +p 2 x , 



(28) 



p x is the electron's longitudinal momentum, the electro- 
static field E x (x) and vector potential a(x) are given by 
Eqs. (9) and (11), respectively, and dotted quantities are 
differentiated with respect to time. 

Equations (26)-(27) can be derived from the Hamilto- 
nian: 



H{x,p x ) = 'f(x,p x ) - (f>{x) , 
where the electrostatic potential reads 



0. 



-,n x* 



x < 0, 

< x < Xb . 



(29) 



(30) 



The Hamiltonian H(x,p x ) is a conserved quantity and 
we can thus write an explicit expression for the electron 
orbit with initial conditions Xo,p x o: 

p x (as) - ± \J[H{x 0l p x0 )+<p{x)f -a?{x)-l . (31) 

Equation (31) suffices to plot portraits of the single- 
electron phase space, as shown in Fig. 2. In the fol- 
lowing subsections we explain how the several solutions 
depicted in Fig. 2 are interrelated, in order to understand 
how phase space geometry affects the threshold of RSIT. 



FIG. 2. (color online) Typical single-electron phase space 
portrait for Eqs. (26)-(27). The first six equilibria 
Qb, Qi, • • • , Qs corresponding to positions Xb, Xi, . . . , £5 and 
zero momentum are shown as blue dots. Separatrices are 
shown as red, dashed lines and some typical trajectories are 
depicted as black, solid lines. The CSL is depicted as a gray- 
shaded area. 



B. Equilibrium solutions 

The simplest type of solutions of Eqs. (26)-(27) are 
equilibrium solutions for which x = p x = 0. We have 
already seen that, within the framework of the stationary 
cold-fluid model, the force balance Eq. (6) is satisfied 
in the plasma and in particular at x = Xb- Thus, the 
point (x,p x ) — (xf,,0) is an equilibrium which we label 
as Qf,. (For the same reasons, any point in the plasma 
with p x — will be an equilibrium.) 

In the CSL and vacuum, on the other hand, the pon- 
deromotive and electrostatic forces are not balanced in 
general, and equilibria for the motion of a test particle 
have to be found by setting x = p x = in Eqs. (26)-(27). 
We label equilibria at the left of Q& as Q m , m = 1, 2, ... , 
where m increases with decreasing x m . 

For x < (in the vacuum) equilibria correspond to 
dxl = -gf = 0, i.e. a — or da/dx = 0, which, according 
to Eq. (11), leads to 



rib 



\/2a 



+ Xb — kir/2 . 



(32) 



Here, k can be any positive integer provided that x k < 0, 
and k even or odd correspond to a{x^) = or a'{x^) = 0, 
respectively. We note that in our labeling scheme, index 
k in xjl does not always correspond to index m in labeling 
of equilibria Q m , i.e. we will generally have x m — x^ with 
rn ^ k. [26] 

For < x < Xb (in the CSL), the equilibrium con- 
dition d x (f> — d x j must be solved numerically, using 
Eqs. (11) and (30) for a(x) and 4>(x), respectively. A per- 
turbative solution can be obtained in the neighborhood 
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of Xb, by expanding p = —g-^ — E x (x) = to second 
order in x — Xb- We obtain two solutions, x = Xb and 



X\ ~ x b + 



2{l + a 2 ) 2 [2{a 2 + n ) - 3n (l + og) 1 ^] 



a?)V2( 4 . 



2a 2 . 



6a 2 ) 



(33) 

Comparing Eq. (33) with condition (20), we see that x\ < 
, as long as a standing wave solution exists, i.e. for 
no > nth- At threshold, uq = n t h, we have x% = xj,. 
That is, if we approach the RSIT threshold (as predicted 
by cold-fluid theory), the equilibrium Qi approaches Q& 
until the two states coalesce, see Fig. 3(c). 



C. Stability of equilibria 

Linear stability analysis of the equilibria determined 
in Sec. Ill B can give us information on the behavior 
of orbits in the neighborhood of the equilibria. For 
notational convenience, we define phase space variables 
C = (Cij C2) = [x-i Px) and rewrite the equations of motion 
[Eq. (26) and Eq. (27)] in the form 

ti=Fi(0- (34) 

where F X (C) = C2/7 and F 2 (C) = -|£ - E x (b). Con- 
sidering infinitesimal perturbations in the neighborhood 
of equilibrium C (m) , and substituting = C (m) + £(t), 
with ||£| <C 1, in Eq. (34), one obtains 

i = A(U)C, (35) 
where the Jacobian matrix A(( m ), with elements 

8R 



(36) 



has been introduced. 

Solutions of the linear system Eq. (35) are of the form 
£(i) = exp[A(( m )t]£,(0), and thus the linear stability of 
equilibrium Q m is determined by the eigenvalues of the 
Jacobian matrix. In Hamiltonian systems with one de- 
gree of freedom, classification of equilibria Q m by linear 
stability is straightforward (see, e.g., Ref. [27]), as there 
are only two possibilities: 

• A(( m ) has a pair of real eigenvalues Ai = — A 2 > 0. 
Solutions then deviate from Q m at an exponen- 
tial rate, ||£(i)|| ~ e Al *||£(0)||, and the equilibrium 
(called a saddle) is unstable. 

• .4 (Cm) has a conjugate pair of purely imaginary 
eigenvalues Xi = = iw. Solutions then oscillate 
around Q m with period 2ir/w, and the equilibrium 
(called a center) is (neutrally) stable. 

Taking into account equilibrium conditions x — p x — 
0, we find from Eq. (36) 



A(Cm) 



l/7„ 
A 2 i 



where 



A21 = 



a 2 (a' ) 2 

u m y^m) 
/ m 



(a'm) 2 + al 



no , x > , 
0, x < 0. 



Here, we have defined a m — a(x m ), a! m = a'(x m ), j m — 
a/1 + a 2 (x m ), and we have used Eq. (8). 
Eigenvalues of A(( m ) are given by 

Ai, 2 (£m) = ±V A 2 \h m . (37) 

In the vacuum, x < 0, equilibria correspond to ei- 
ther a(x^) = (k even, nodes of the standing wave) 
or a'(x7) — (k odd, antinodes of the standing wave), 
where the x^ are given by Eq. (32). Then, Eq. (37) yields 
by using Eqs. (32) and (11), 



(iy/2a , 

I V 1+2a o 



k even, 

k odd. (38) 



Thus, in the vacuum, equilibria alternate between being 
(neutrally) stable (k even, nodes) and unstable (k odd, 
antinodes). 

In the CSL, x > 0, we have 



Ai, 2 (^m) = i-FV^mK™ - ( a 'm) 2 ) + a 2 m {a' m ) 2 - t| 



In 



no ■ 



For the equilibrium Qh at the plasma boundary, x — 
Xb, we get from Eq. (11) 

Ai |2 (i 6 ) =±Ja% + 2a 2 - 2a 2 - n Q (l + a 2 )V 2 /{l + a 2 b ) . 

(39) 

Linear (neutral) stability of Q b requires 

at + 2a\ - 2al - n (l + a 2 b f /2 < , 

or, using Eq. (18) to eliminate ao, 

2(a 2 b + n ) ~ 3n Q (l + a 2 b )^ 2 <0. (40) 

The same condition for linear stability of the equilibrium 
at Xb was obtained by Eremin et al. [19] by considering 
the infinitesimal variation in electrostatic and pondero- 
motive force experienced by an electron whose position 
has been perturbed infinitesimally to Xf,— \8x\. Condition 
(40) also coincides with condition (20) of existence of a 
stationary standing wave obtained by Cattani et al. [17]. 
Therefore, as long as an equilibrium at x b exists, it is 
neutrally stable. 

Assessing stability of the equilibria with < x\ < x b 
analytically is somewhat more difficult [even when an 
explicit expression such as Eq. (33) is available]. We 
can, however, conclude that Qi is an unstable equilib- 
rium on topological grounds. If we assume Qi to be 
stable, then motion in its neighborhood would be oscil- 
latory. Therefore, a point (x s ,0) in phase space with 
x\ < x s < x b would be shared by oscillatory solutions 
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FIG. 3. (color online) For a given laser field amplitude, here ao = 7, the absolute value of the critical momentum for an 
electron to escape to the vacuum \p"\ decreases as no decreases. Shown are the cases (a) no = 5.8, (b) no = 4.8, (c) no = 3.31. 
Color-code is the same as described in the caption of Fig. 2. In panel (c) equilibria Qi and G> cannot be distinguished within 
the resolution of this plot, no being slightly above the cold-fluid theory threshold nth = 3.30458. Note that the scale of x and 
p x has been kept the same in all panels. 



encircling Qi and Qb (in phase space). This would con- 
tradict uniqueness of solutions, unless the point (x s ,0) 
were to be reached in infinite time, i.e. unless it is an un- 
stable equilibrium. However, by construction there is no 
equilibrium beween Q\ and Qb- In fact, the degenerate 
oscillations introduced in this argument, which reach Q\ 
in infinite time, are the familiar separatrices of bounded 
and unbounded motion, which we will now study in de- 
tail. 



D. Separatrices 

In the vacuum (x < 0), all unstable equilibria at x^ 
(with k odd) correspond to the same value of H, 



2a 2 . 



(41) 



Conservation of H, thus allows for a heteroclinic connec- 
tion, i.e., for an orbit which starts infinitesimally close 
to Qfe and ends infinitesimally close to Qk+2 or Qfc-2 (in 
infinite time). According to Eq. (31) these orbits obey 



p x (x) = ±yj: 



2a\ 



a 2 (x) 



(42) 



Heteroclinic connections, Eq. (42), act as separatrices of 
bounded and unbounded motion, see Fig. 2. 

Within the CSL (0 < x < Xb), an unstable equilibrium, 
e.g. Qi in Fig. 2, will in general have H(xi,0) ^ H{xz, 0) 
since H now also includes an electrostatic field contribu- 
tion. Therefore, a heteroclinic connection from Q3 to Qi 
is not possible, and the separatrix starting out at Q3 is a 
homoclinic connection, i.e. an orbit that returns to Q3 in 
infinite time. For the same reason, the separatrix labeled 
B in Fig. 2 starts in the neighborhood of Qi and wanders 
off to x = —00, while the separatrix labeled A starts at 
x = —00 and ends at Qi. 

Of greatest importance in the following discussion are 
the separatrices labeled T and A, as they determine the 



region within which motion close to Qb is oscillatory. 
The equations of the separatrices T and A are given by 
Eq. (31) with (xo,p x o) — (xi, 0) [on separatrix T, motion 
is backwards in time and (xi,0) is a final, rather than 
initial, condition]. The point on separatrix T at posi- 
tion Xb (at the plasma boundary) then defines a critical 
momentum p", given by 

-1 1/2 



^l + a 2 (x 1 ) + n Q (x 2 1 -xl)/2 - a\ - 1 

'(43) 

If a single electron at the edge of the plasma Xb is 
given an initial momentum — |Ap x |, with |Apa;| < \p x r \, 
it will move within the limits set by separatrices T and 
A, returning back to the plasma. If, on the other hand 
|Ap x | > \p"\, the electron's motion will be unbounded 
and it will escape to the vacuum. Alternatively, one can 
define a critical value of the Hamiltonian 



H CT = H {x x , 0) = v/l + a 2 (xi) + n x\/2 . 



(44) 



Motion of electrons with H(xb,p x ) > H cl ' and p x < 
will be unbounded. 

Equation (43) shows that \p"\ is always non-zero as 
long as xi ^ Xb; for fixed ao it becomes smaller as no 
decreases and x\ approaches Xb, vanishing at the thresh- 
old n t h given by Eq. (20). This behavior is illustrated in 
Fig. 3 for a = 7. (See also Fig. 10.) 

With the above results it becomes clear that finite per- 
turbations of initial conditions of electrons at the edge 
of the plasma, for example due to longitudinal electron 
heating, could lead to electrons escaping toward the vac- 
uum even when Qb is stable in the linear approximation, 
provided that the perturbation (here negative momen- 
tum) is large enough. Our main conclusion is that pulse 
propagation by expulsion of electrons toward the vacuum 
could occur for densities higher than the threshold den- 
sity nth predicted by the cold fluid approximation. In 
Sec. IV we show that electron heating at the edge of the 
plasma indeed provides a mechanism by which electrons 



acquire sufficient momentum to escape toward the vac- 
uum. 

IV. PIC SIMULATIONS 

To investigate the transition from total reflection 
to RS1T, we perform PIC simulations [28] using the 
one-dimensional in space, three-dimensional in velocity 
(1D3V) code Squash [29]. The code uses the finite- 
difference, time-domain approach for solving Maxwell's 
equations [30] , and the standard (Boris) leap-frog scheme 
for solving the macro-particle equations of motion [31]. 
Charge conservation is ensured by using the method pro- 
posed by Esirkepov when projecting the currents [32]. 

In all simulations presented here, ions are immobile 
and only electron motion is considered. We use the spa- 
tial resolution dx = Al/500 and time step dt = rr,/1000, 
where Al and tl are the laser wavelength and duration 
of one optical cycle, respectively. Up to 1000 macro- 
particles per cell have been used. 

The plasma extends from x — to x — L p , with a 
constant initial density no and electron temperature Tq ~ 
5T0~ 4 (in units of m e c 2 ). The plasma size L p is chosen so 
that L p > c Tint , where Tint is the laser-plasma interaction 
time. Hence, the plasma is long enough to be considered 
semi-infinite. The CP laser pulse [as described by Eq. (3)] 
is incident from x < onto the plasma. In this work we 
consider laser field amplitudes in the range ao — 1 — 30. 
The laser pulse profile is trapezoidal, i.e. the intensity 
increases linearly within a rise time T r , up to a maximum 
value clq/2, and we consider the exemplary cases r r — 
0.25 and t t — At^. 

Figure 1(a) summarizes our findings on RSIT, compar- 
ing the threshold density nth( a o) predicted by Cattani et 
al. [17] with our 1D3V PIC simulation results. In order 
to determine whether RSIT occurs or not in a simula- 
tion, the position Xb of the maximum electrostatic field 
is plotted as a function of time (see Fig. 4). The regime 
of total reflection is characterized by the formation of a 
CSL with (approximately) constant thickness Xb (Fig. 4, 
no = 6.7—8). On the other hand, RSIT is associated with 
front propagation at an approximately constant velocity 
Vf, so that the position of the maximum electrostatic 
field increases linearly with time (Fig. 4, no — 5.75 — 6). 
This allows us to place lower and upper bounds on RSIT 
threshold density, for a certain a , indicated by error bars 
in Fig. 1. For densities within these limits, it is hard to 
decide whether RSIT occurs or not (Fig. 4, n = 6.25). 

In the next subsections we examine in detail typical 
cases of total reflection and front penetration. 



A. Total reflection 

Whenever total reflection occurs, the system eventu- 
ally settles to a quasi-stationary state. The size of the 
charge separation layer Xb remains constant or slightly 




FIG. 4. (color online) Position of the maximum value of the 
electrostatic field x\, as a function of time from PIC simula- 
tions with different densities and ao = 15, t t = 0.25 tl- 



oscillatory around a value that is found to be in good 
agreement with the theoretical prediction of the cold- 
fluid model [Eq. (17)], see Fig. 5. The same is true for 
the field and density profiles; a worst case agreement is 
shown in Fig. 6, where the quasistationary state reached 
for a — 15, no — 1 and r r = 0.25r^ is close to the nu- 
merical RSIT threshold (the agreement becomes better 
for higher no or larger r r ). Although the density profile 
presents oscillations, the fields in the CSL and vacuum 
agree very well with the predictions of cold-fluid theory. 
This justifies a posteriory our use of stationary cold fluid 
theory predictions for the fields in the vacuum to ana- 
lyze single electron phase space in Sec. III. The phase 
portrait for ao = 15, no = 7 and r r = 0.25tx is shown 
in the top row of Fig. 7. It is clearly seen that electrons 
in the CEL do not have zero longitudinal momentum p x 
as the stationary cold-fluid model suggests, but rather 
oscillate around Xb [the latter being in good agreement 
with Eq. (17)]. As the minimum momentum attained by 
electrons, which we will call p™ m , is smaller in absolute 
value than the critical momentum required to move be- 
yond the limits set by the separatrices of bounded and 
unbounded motion, |p™ ln | < \p x r \, electrons which cross 
the plasma boundary Xb do not escape into the vacuum 
but rather re-enter the CEL. 

As can be seen in Fig. 4, for no = 6.7 — 8, the position 
of the plasma boundary Xb oscillates in time, leading to 
oscillations of the maximum electrostatic field. These os- 
cillations can be related to the excursion of electrons in 
the region x < Xb, cf. the top panel of Figure 7. To verify 
this, we plot in Fig. 8 the period T osc of these oscillations 
for different ao and n well in the regime of total reflec- 
tion. The frequency of these oscillations is not linked to 
the plasma frequency (observe the dependence on ao in 
Fig. 8) but rather on the frequency of oscillations of elec- 
trons around the equilibrium Qb- If we ignore the role 
of the self-consistent fields within the plasma, the char- 
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FIG. 5. (color online) Comparison of cold-fluid model predic- 
tion for Xb (blue, solid lines) with the (time-averaged) posi- 
tion of the maximum electrostatic field in our PIC simulations 
(dots), with r r = 0.25 tl. For values of no to the left of the 
thick, gray, dashed line RSIT occurs and Xb does not reach a 
constant average value in our PIC simulations. 
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FIG. 6. (color online) Electron density and field profiles from 
PIC simulations for ao = 15, no = 7 at t — 2.55 tl (top 
panel) and t — 2.95 tl (bottom panel). The stationary cold- 
fluid model solution for the electron density (red, dashed line), 
electrostatic field (blue, dotted line) and vector potential en- 
velope (black, dash-dotted line) are also shown. Note that 
densities have been rescaled to the unperturbed density no 
for better readability. 



acteristic period of oscillation in the linear neighborhood 
of Qb reads Tq 6 = 27r/ImAi, where Ai is the eigenvalue 
given by Eq. (39) . As shown in Fig. 8, we find T osc cx TQ b . 
We also note the similarity of these oscillations with the 
so-called piston oscillations in laser hole-boring [10], al- 
though in the present case the oscillations only involve 
electrons. 



B. RSIT 

The cold fluid model presented in Sec. II predicts a 
sharp threshold, either for density uq or laser amplitude 
ao, for RSIT. However, as already mentioned above, one 
of the main results of this paper is that our PIC simu- 
lations clearly show RSIT in a parameter region where 
the cold fluid model predicts total reflection [area (B) 
in Fig. 1(a)]. A typical case of RSIT in this regime 
is presented in Fig. 9, where ao = 15, uq = 5.5 and 
r r = 0.25t£. Charge separation and compressed elec- 
tron layers are formed in the early stages of interaction, 
with profiles that agree well with the predictions of cold- 
fluid theory. However, electrons escape the CEL, and 
the pulse can propagate (see middle row of Fig. 9). The 
mechanism of propagation is rather complex, but its ini- 
tial phase can be intuitively understood as follows. When 
a sufficiently high number of electrons escapes from the 
CEL to the vacuum, the electrostatic field within the CSL 
decreases, the ponderomotive force is no longer balanced 
and the laser pulse can push the CEL deeper into the 
plasma. The increase of the CSL size tends to compen- 
sate the force imbalance, but as more and more electrons 
escape, the pulse continues to propagate deeper into the 
plasma. We note that once electrons escape and propaga- 
tion commences the stationary model is no longer valid 
and electron dynamics becomes complex, with electron 
bunches leaving and re-entering the plasma (see Fig. 9 
and Ref. [19]). 

To understand how the shrinking of the width of sepa- 
ratrices in phase space with decreasing density (and con- 
stant ao) leads to propagation, we examine the phase 
space portrait for ao = 15, no — 6.0 and r r = 0.25 Tl, 
which corresponds to a case just below the numerical 
density threshold for RSIT, see the bottom row of Fig. 7. 
In this case, the minimum momentum acquired by elec- 
trons in the CEL satisfies \p™ m \ > \p"\ and electrons 
move outside the separatrix of bounded and unbounded 
motion, eventually reaching the vacuum, while the CEL 
moves deeper into the plasma. 

Figure 10 provides a further verification of the role the 
longitudinal electron heating plays in enabling electrons 
to escape from the CEL into the vacuum. We use Eq. (43) 
to plot |p"'| as a function of ao and n (light-gray sur- 
face). For a given rise time, here T r = 0.25 tl, we also 
plot, as a function of ao and no, the absolute value of 
the minimum momentum |p™ m | acquired by electrons in 
the CEL as inferred from our PIC simulations (dark-blue 
surface). To reduce noise we average |p™ m | over one laser 
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FIG. 7. (color online) Comparison between phase space separatrices as predicted by the stationary, cold-fluid model, and single 
particle distribution function f(x,p x ) from PIC simulation results for ao = 15 and no = 7 (top row), no = 6 (bottom row) and 
rise time r r = 0.25 tl- Snapshots are shown 0.2 tl apart. The plasma boundary (x — Xb), as predicted by the cold-fluid model, 
is indicated by a black, dotted, vertical line. The color coding of trajectories follows Fig. 2. Note the logarithmic scale in the 
color coding of f(x,p x ). 
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FIG. 8. (color online) Period of oscillations of xt as a function 
of ao for three different no. Points with error bars correspond 
to the periods as deduced from our PIC simulations, while the 
solid lines correspond to T osc = 1.75 Tq 6 = 1.75 x 27r/ImAi, 
where Ai is given by Eq. (39). 



period (starting at t ~ 2rr,), or at most until electrons 
escape. Thus, our |p" lln | is generally slightly underesti- 
mated, however the intersection of the two surfaces \p"\ 
and |p™ m | lies within the limits set by the error bars in 
Fig. 1(a). Note that |p™ ln | is getting smaller with decreas- 
ing ao, and one recovers the threshold predicted by cold- 
fluid theory for a < 5, where the longitudinal electron 
momenta become negligible [compare with Fig. 1(a)]. 



C. Dependence on rise time 

As we have seen, the threshold for transition between 
total reflection and RSIT clearly depends on the longi- 
tudinal momenta of the electrons in the CEL. As these 
momenta come from collisionless heating of the electrons, 
we may expect that the RSIT threshold also depends on 
the laser pulse profile. As can be seen in Fig. 1(a), the de- 
viation of the numerically obtained RSIT threshold from 
the predictions of cold fluid theory is smaller for a pulse 
with larger rise time, suggesting a less significant electron 
heating in the CEL at given ao and uq. The effect of pulse 
rise time on the width of the longitudinal electron mo- 
mentum distribution function is shown in Fig. 11, where 
the space-integrated distribution for ao = 15, no = 7 is 
compared for the cases r r = 0.25t l and r r = Ar L . The 
stiffcr pulse clearly corresponds to a larger |p™ m |. 

In Fig. 12, we moreover compare the front propaga- 
tion velocity, Vf, for two sets of simulations with rise 
times tv = 0.25tx and t t — Att,- The front propaga- 
tion speed Vf is determined by the slope of the curves 
Xb(t), see Fig. 4. We have studied cases of propagation 
for different ao and no close to the threshold predicted 
by cold- fluid theory, for which Vf ranges from 10~ 3 c up 
to 0.25 c, see Fig. 12. Within the error bars for the trans- 
parency threshold, Vt takes values too small to reliably 
indicate propagation (i.e. beyond the accuracy permit- 
ted by our spatial and temporal resolution). As Fig. 12 
shows, the propagation velocity v/ for the same ao and 
no is generally lower for the pulse with the larger rise 
time, r r = 4tx. Nevertheless, for higher ao, Vf is far 



11 




E x (x) 


n e (x) 


m(x) — -^pr 


" 










10 -5 



10 -5 



10 



FIG. 9. (color online) Field and density evolution for RSIT above the cold-fluid theory threshold, ao = 15, no = 5.5. Snapshots 
are taken one laser period apart, starting at t = 0. 



from negligible for densities lying well above the cold- 
fluid threshold, even for the case with larger rise time. 



V. DISCUSSION AND CONCLUSIONS 

The relativistic, cold-fluid, stationary solutions of 
Refs. [15, 17, 18] provide a convenient starting point to 
investigate the threshold of RSIT, even in the presence 
of longitudinal electron heating. While the fields inside 
the plasma clearly differ from the predictions of cold fluid 
theory, the fields in the CSL and vacuum are rather insen- 
sitive to density fluctuations within the plasma. There- 
fore, the dynamics of a test electron in the CSL or the 
vacuum can be accurately described using the fields of 
the stationary problem. This finding allows us to specify 
separatrices of bounded and unbounded motion for sin- 
gle electron dynamics, encapsulating the competition of 
ponderomotive and electrostatic forces at the edge of the 
plasma. 

We have shown that one can define a critical momen- 
tum \p™\, Eq. (43), or value of the Hamiltonian H cl , 
Eq. (44), corresponding to the separatrix which delim- 
its oscillatory motion around the equilibrium position at 
the edge of the CEL. When a sufficiently high number of 



electrons at the edge of the CEL have \p x \ > \p"\ and 
escape to the vacuum, RSIT occurs. In this work, we 
did not focus on the mechanism that provides momen- 
tum to electrons, i.e. we did not attempt to provide a 
model for the collisionless heating mechanism. We did 
however show, through our numerical study of the im- 
pact of the pulse rise time, that the pulse shape crucially 
affects longitudinal heating and that stronger heating re- 
sults in a higher threshold density for RSIT. A detailed 
model for electron heating, which would allow us to pre- 
dict \p™ in \ rather than infer it from PIC simulations, as 
done in Fig. 10, will be pursued elsewhere. We stress 
that, although in more realistic scenarios of laser-plasma 
interaction the actual heating at the plasma boundary 
would depend on several factors (see Ref. [33] for a re- 
cent study), the basic mechanism of electron escape into 
the vacuum at high enough momentum is expected to be 
the same. 

In summary, we have used a dynamical systems ap- 
proach to bridge the cold-fluid and kinetic levels of RSIT 
description. Deviations of PIC simulations from cold- 
fluid theory predictions are explained as a longitudinal 
heating effect induced by the incident laser pulse. The 
pulse temporal profile clearly affects electron heating and 
through it the threshold of RSIT. While there are sev- 
eral experimental works addressing RSIT in the case of 
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FIG. 10. (color online) The absolute value of momentum \p c J\ 
corresponding to the separatrix of bounded and unbounded 
motion for electrons at xi for different ao and no, according to 
Eq. (43) is shown as a light-gray surface. An estimate of the 
absolute value of the minimum momentum |p™ m | attained by 
electrons in the CEL, as determined by our PIC simulations 
with r r = 0.25tl, is shown as a dark-blue surface. The light- 
arid dark-colored points represent PIC simulation results cor- 
responding to onset of RSIT and total reflection, respectively, 
for ao = 5, 10, . . . , 30. On both surfaces, lines of constant ao 
are drawn to guide the eye. The contour \p"\ = (black, 
thick, solid line) corresponds to the threshold for RSIT pre- 
dicted by cold-fluid theory. 




FIG. 12. (color online) Front velocity t>/ as measured from 
PIC simulations with two different pulse rise times, (a) tv = 
4tx, (b) 7y = 0.25 tl- The blue, solid line and error bars 
are the same as described in the caption of Fig. 1(a). The 
lower density range of simulations performed for a given ao 
has been set to improve readability. 




linearly polarized laser pulses [34-37], to the best of our 
knowledge the verification of RSIT for CP light remains 
elusive. We hope that our results trigger further inves- 
tigations in this domain, as the reported dependency of 
the RSIT threshold on the pulse profile could provide a 
versatile tool for high-contrast CP laser pulse character- 
ization. 



FIG. 11. (color online) Space-integrated (over all x) longitu- 
dinal momentum distribution f(p x ) at t = 15 tl for ao = 15, 
no = 7 and r r — 0.25 tl (red, solid curve), r r — 4tl (black, 
dotted curve). Both cases correspond to total reflection. The 
central peak corresponds to the plasma bulk. 
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